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SINGULAR  INTEGRALS,  IMAGE  SMOOTHNESS,  AND  THE 
RECOVERY  OF  TEXTURE  IN  IMAGE  DEBLURRING 

ALFRED  S.  CARASSO* 


Abstract.  Total  variation  (TV)  image  deblurring  is  a PDE-based  technique  that  preserves 
edges,  but  often  eliminates  vital  small-scale  information,  or  texture.  This  phenomenon  reflects  the 
fact  that  most  natural  images  are  not  of  bounded  variation.  The  present  paper  reconsiders  the  image 
deblurring  problem  in  Lipschitz  (Besov)  spaces  A (a,p,q),  wherein  a wide  class  of  non-smooth  images 
can  be  accomodated,  and  develops  a fast  FFT-based  deblurring  method  that  can  recover  texture  in 
cases  where  TV  deblurring  fails  completely.  Singular  integral  modifiers,  such  as  the  Poisson  kernel, 
are  used  to  create  an  effective  new  image  analysis  tool  that  can  calibrate  the  lack  of  smoothness  in  an 
image.  It  is  found  that  a rich  class  of  images  belong  to  A(a,  1,  oo)D  A(/3,  2,  oo),  with  0.2  < a,  /3  < 0.6. 
The  Poisson  kernel  is  then  used  to  regularize  the  deblurring  problem  by  appropriately  constraining 
its  solutions  in  A(a,2,  oo)  spaces,  leading  to  new  L2  error  bounds  that  substantially  improve  on  the 
Tikhonov-Miller  method.  This  new  so-called  Poisson  Singular  Integral  or  PSI  method  is  found  to 
be  well-behaved  in  both  the  L 1 and  L2  norms,  producing  results  closely  matching  those  obtained 
in  the  theoretically  optimal,  but  practically  unrealizable  case  of  true  Wiener  filtering.  Deblurring 
experiments  on  synthetically  defocused  images  illustrate  the  PSI  method’s  significant  improvements 
over  both  the  total  variation  and  Tikhonov-Miller  methods. 

Key  words,  image  deblurring,  singular  integrals,  non-smooth  images,  total  variation,  loss 
of  texture,  Lipschitz  spaces,  Besov  spaces,  Poisson  kernel,  semi-group  approximation,  recovery  of 
texture,  Tikhonov-Miller  method,  true  Wiener  filtering,  PSI  method. 
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1.  Introduction.  The  space  BV(R 2 ) of  functions  of  bounded  variation,  normed 
by  the  ‘total  variation’  seminorm  jR2  | V f\dxdy,  plays  an  important  role  in  much  recent 
work  in  image  analysis.  See  e.g.,  [11],  [13],  [14],  [15],  [16],  [20],  [21],  [26],  [33],  and  [39]. 
In  particular,  highly  successful  applications  of  the  total  variation  approach  to  image 
denoising  have  been  well-documented.  In  contrast,  total  variation  image  deblurring  is 
generally  not  well-behaved,  and  often  results  in  unacceptable  loss  of  fine  scale  infor- 
mation. This  phenomenon  is  now  believed  traceable  to  an  improper  choice  of  function 
space  [24],  The  present  paper  reconsiders  the  image  deblurring  problem  in  Lipschitz 
spaces  A (a,p,q)  wherein  a wide  class  of  non-smooth  images  can  be  accomodated. 
A new  and  fast  FFT-based  deblurring  technique  is  developed  that  can  demonstrably 
recover  texture  in  cases  where  total  variation  deblurring  fails  completely.  The  approx- 
imation properties  of  certain  singular  integral  modifiers  are  intimately  linked  to  such 
Lipschitz  spaces  [3],  [4],  [36].  Here,  these  properties  are  exploited  in  two  distinct  ways. 
In  the  first  part  of  the  paper,  such  singular  kernels  are  used  to  create  an  effective  new 
FFT-based  image  analysis  tool  that  can  calibrate  the  lack  of  smoothness  in  an  image. 
This  tool  can  be  used  in  contexts  unrelated  to  deblurring,  e.g.,  as  a sharpness  analysis 
tool  in  performance  evaluation  of  imaging  systems  or  image  reconstruction  software 
[38],  or  as  a tool  for  detecting  and  quantifying  ‘fine-structure’  content  in  images.  In 
the  second  part  of  the  paper,  singular  integral  modifiers  are  used  as  regularization 
tools  in  the  deblurring  problem.  Specifically,  we  show  how  to  stabilize  ill-posedness 
by  using  the  Poisson  kernel  to  impose  a-priori  constraints,  in  appropriate  A (a,p,oo) 
spaces,  on  the  desired  non-smooth  deblurred  image.  This  so-called  Poisson  Singular 
Integral  or  PSI  method,  is  only  one  of  a wide  variety  of  singular  integral  deblurring 
methods  that  can  be  constructed.  Restricting  attention  to  the  case  of  defocus  blurs, 
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we  derive  L1  error  bounds  for  the  PSI  method  for  images  € A(q,  2,  oo),  and  demon- 
strate robust  recovery  of  fine  structure  in  synthetically  blurred  images.  Wider  classes 
of  blurs  will  be  considered  in  future  papers. 

Extensive  numerical  experiments  with  known  exact  solutions  indicate  that  the 
PSI  method  is  remarkably  well-behaved.  In  both  the  Ll  and  L2  norms,  relative 
errors  in  the  PSI  method  are  found  to  closely  approximate  those  obtained  in  the 
theoretically  optimal,  but  practically  unrealizable  case  of  true  Wiener  filtering.  The 
latter  method  requires  prior  knowledge  of  the  exact  power  spectra  of  each  of  the  noise 
and  the  unknown  desired  sharp  image,  i.e.,  a total  of  8N'2  prior  data  values  for  a 
2N  x 21V  image.  Since  the  PSI  method  requires  only  4 prior  data  values,  its  ability 
to  closely  track  Wiener  filtering  is  especially  noteworthy.  The  availability  of  reliable 
fast  deblurring  methods  is  of  major  significance  in  the  case  of  non-smooth  images. 
The  true  value  of  the  Lipschitz  exponent  a in  the  desired  sharp  image  is  usually  not 
known  in  advance,  although  a plausible  range  of  values  for  a can  often  be  deduced. 
Fast  algorithms  enable  simultaneous  computation  and  display  of  large  numbers  of  trial 
deblurred  images,  resulting  from  multiple  choices  for  a,  or  some  of  the  regularization 
parameters.  We  stress  that  the  PSI  method  is  exclusively  intended  for  deblurring  and 
is  not.  intended  for  denoising. 

2.  Lack  of  smoothness  of  images.  In  [27],  a new  analytical  framework  for 
image  processing  is  introduced,  in  which  a given  image  f(x,y)  is  conceptualized  as 
being  the  sum  of  three  components,  f(x,y)  = u(x,y)  T-  v(x,y)  + w(x,y).  Loosely 
speaking,  u(x,y)  contains  the  edges  and  the  other  high-priority  information  that  is 
sufficient  for  object  recognition,  v(x,  y)  contains  the  fine-scale  details  and  other  low- 
priority  information  that  is  often  not  necessary  for  recognition,  and  w(x,  y)  represents 
noise.  The  v(x,y)  component  is  called  texture.  One  example  of  v(x,y)  might  be  the 
hair  in  a photograph  of  a person’s  face.  Another  example  of  v(x,  y)  might  be  the  heat- 
shield  tiles  in  an  image  of  the  Columbia  space  shuttle.  The  ability  to  resolve  individual 
hairs  is  generally  not  necessary  for  identification.  In  several  image  processing  tasks, 
such  as  compression,  segmentation,  or  face  recognition,  this  texture  component  can 
often  be  neglected.  However,  there  are  numerous  other  situations  where  v(x,y)  may 
be  of  paramount  interest.  It  is  shown  in  [27]  that  only  the  u(x,  y)  component  can 
generally  be  expected  to  lie  in  BV(R2).  In  [24],  it  is  proved  that  most  natural  images 
are  not  of  bounded  variation,  because  the  texture  component  v(x,y)  generally  has 
infinite  total  variation. 

Denoising  and  deblurring  are  two  basic  image  processing  tasks  where  total  varia- 
tion restoration  has  been  extensively  applied.  Such  restoration  can  be  accomplished 
most  effectively  by  solving  an  initial  value  problem  for  an  appropriate  nonlinear 
anisotropic  diffusion  equation,  using  the  stepwise  marching  scheme  described  in  [26]. 
In  deblurring,  one  typically  starts  with  a degraded  image  g(x,y)  which  differs  from 
the  desired  true  image  f(x,y)  in  that  the  u(x,y)  component  is  blurred  but  recogniz- 
able, the  v(x,  y)  component  is  seriously  attenuated  and  often  not  recognizable,  and 
the  w(x,y)  component  is  usually  small.  Reconstructing  v(x,y)  while  keeping  w(x,y) 
small,  is  the  prime  objective  in  numerous  medical,  astronomical,  industrial,  and  sci- 
entific contexts  [9],  [10].  However,  while  total  variation  deblurring  sharpens  u(x,y) 
and  keeps  w(x,y)  small,  the  texture  component  v(x,y)  is  often  eliminated  due  to  the 
‘staircase  effect.’  [13],  [20],  [29],  [30],  [39].  This  is  in  accordance  with  the  analyses  in 
[24],  [27], 

Let  x = (, Xi,X2 ) € R2.  Postulating  f(x)  € BV(R 2)  means  that  f(x)  is  con- 
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strained  to  satisfy, 


(1)  / \f  (x  + h)  — f(x)\dx  < Const  \h\. 

J R2 

However,  from  the  standpoint  of  modeling  texture,  it  is  advantageous  to  consider 
functions  f(x)  satisfying  weaker  constraints,  such  as 


/r2 


| f(x  + h)  — f(x)\pdx  > < Const  |h| 


i/p 


0 < a < 1. 


Such  an  / lies  in  A(a,p,oo).  With  0 < a < 1,  1 < p < oo,  the  Lipschitz  (Besov) 
spaces  A (a,p,q)  [36],  [37],  consist  of  the  class  of  functions  f(x)  € LP(R 2)  with  finite 
seminorm  ||  / ||Qpg,  where 


(3)  ||  / 


\apq 


f(x  + h)  — f (x)  \\py  dh/\h\2 


1/9 


1 < q < oo, 


(4) 


/ lUpoo  = sup  {\h\  a ||  f(x  + h)  - f(x)  ||p}  , q = oo. 
h€R2 


For  given  p and  q , functions  with  larger  values  of  a are  better  behaved,  or  ‘smoother’, 
than  functions  with  smaller  values  of  a,  and  functions  in  A(a,p,q\)  are  smoother  than 
functions  in  A(a,p,  <72)  if  qi  < <72  • In  fact,  the  following  continuous  embedding  results 
are  proved  in  [37,  Theorem  9|. 


(5)  A(q2,p,  Qi)  C A(qi,p,  (72)  0 < cti  < Q'2  < 1;  1 < qi  < qo  < 00. 

Also,  in  R2, 


(6)  A(a,p,q)  C A((3,r,q),  a - 2/p  = (3 -2/r,  p < r. 

Let  r = 2,  let  the  pair  (ct,p)  satisfy*  2/(1  + a)  < p < 2,  and  let  (3  = 1 + a — 2/p.  Then, 
0 < P < a,  and  it  follows  from  (5)  and  (6)  that 

(7)  A (a,p,  q)  C A(/3, 2,  q)  C A(/3, 2,  00)  C L2(R2). 


This  result  will  be  important  in  the  sequel. 

For  given  fixed  p with  1 < p < 00  and  q — 2 p/  (2  + ap),  a class  of  Lipschitz  spaces 
A(q,  q,  q)  C LP(R2)  is  considered  in  [18],  [19],  and  shown  to  contain  common  types 
of  images.  A method  for  empirically  estimating  image  smoothness  is  developed  in 
[18],  [19],  based  on  analyzing  the  behavior  of  lossy  wavelet  compression  of  the  image 
f(x,y).  In  [12],  the  spaces  A (a,q,q)  C L2(R2),  q = 2/(1  + ct),  are  advocated  as 
being  particularly  appropriate  for  accomodating  a rich  variety  of  real  images  in  an 
L2  setting.  Lossy  wavelet  compression  is  again  used  to  estimate  image  smoothness, 
and  values  of  a in  the  range  0.4  < a < 0.75  are  reported  in  [12]  for  a class  of  24 
test  images  e A(a,  y~,  jfyy).  In  this  approach,  the  restriction  on  q does  not  allow 
consideration  of  the  wider  spaces  A(a,p,oo)  D A (a,p,q).  Clearly,  such  a values  are 
an  indication  of  true  image  smoothness  only  when  the  image  is  largely  noise  free.  If 
the  noise  component  w(x,  y)  is  not  sufficiently  small,  artificially  low  values  of  a must 
be  expected. 
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The  present  independent  method  of  estimating  image  smoothness  rests  on  an  en- 
tirely different  analytical  basis,  and  requires  neither  wavelet  expansions  nor  image 
compression.  Instead,  the  method  uses  fast  FFT  algorithms  to  convolve  the  image 
with  a specific  type  of  kernel,  and  then  analyzes  how  well  this  convok  ed  image  ap- 
proximates the  original  image  as  the  kernel  approaches  the  Dirac  S— function.  This 
simple  direct  approach  permits  consideration  of  the  spaces  A(a,p,  oo),  1 < p < oo. 
The  results  obtained  here  are  compatible  with  those  obtained  in  [18],  [19],  [12],  and 
[24] . We  indeed  find  that  most  natural  images  are  not  of  bounded  variation,  and  that 
a wide  variety  of  images  € A(a,  1,  oo)  with  0.2  < a < 0.6. 

Remark  1.  We  deal  with  high  resolution  images  f(x,  y)  of  size  512x512  or  1024x  1024 
pixels.  Such  an  f(x,  y ) may  be  viewed  as  a piecewise  constant  or  trigonometric  poly- 
nomial approximation  to  the  original  intensity  field  f°°(x,  y),  or  as  some  other  kind  of 
finite  dimensional  representation  of  the  infinite  dimensional  object  f°°.  All  norms  are 
equivalent  on  a finite  dimensional  space.  Hence,  even  if  f°°(x,y ) is  not  of  bounded 
variation,  the  discrete  total  variation  norm  for  f(x,  y)  is  always  finite,  though  it  may 
be  very  large.  To  estimate  smoothness  properties  of  f°°(x,y)  by  examination  of  the 
finite  dimensional  representation  f(x,y)  requires  some  sagacity.  In  [18,  §4B,  §5B],  the 
authors  stress  that  in  their  method  of  estimating  the  value  of  a by  monitoring  the 
rate  of  convergence  as  a function  of  the  number  J\f  of  nonzero  wavelet  coefficients,  one 
must  restrict  attention  to  low  values  of  J\f.  At  high  values  of  J\f,  the  fact  that  f(x,y) 
is  actually  piecewise  constant  causes  the  error  to  decrease  much  too  rapidly,  resulting 
in  an  artificially  high  reading  for  a that  diverges  from  true  behavior  in  f°°(x,  y).  This 
same  finite  dimensionality  pitfall  occurs  in  the  present  approach,  but  wears  a different 
guise.  See  Remark  2 and  the  discussion  surrounding  Figure  1 in  Section  5 below. 

We  shall  use  the  spaces  A(a,  l,oo)  and  A(a,  2,  oo)  for  examining  and  classifying 
image  smoothness.  However,  deblurring  applications  will  be  limited  to  the  spaces 
A(/3,  2,  oo)  C L2(R2),  wherein  all  spaces  A(a,p,  q),  2/(1  + a)  < p < 2,  1 < q < oo, 
are  continuously  embedded.  The  spaces  A(a,  2,  oo)  will  be  shown  to  contain  a rich 
and  significant  class  of  images. 

3.  The  spaces  A (a,p,q)  and  the  Poisson  singular  integral.  Define  the 
Fourier  transform  h(f,r/)  of  h(x,y)  € Z, 1 ( ) bv 

(8)  iF{h}  — h(f,r])  = ( h(x,  y)e~2m^x+vy^  dxdy . 

Jr 2 

For  each  fixed  t > 0,  consider  the  Poisson  kernel  in  R 2 

(9)  i’(x,y,t)  = ,ov3/o.  (x,y)GR2. 

2tt(x-  + y-  + 

We  have 

(10)  . =e~tp,  p=  yjf2  + rf2. 

For  each  t > 0,  define  the  linear  operator  Ut  on  LP(R2),  1 < p < oo,  by 

(11)  Utf=  / ip(x,y,t)f(x  - u,y  - v)dudv 

Jr 2 

It  can  be  shown  that  linpjo  ||  Ul f — / ||p=  0.  Defining  U°  to  be  the  identity  op- 
erator, it  follows  that  for  s,t  > 0,  L rtUs  = Ut+S.  In  fact,  is  a holomorphic 
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contraction  semigroup  on  LP(R2).  See  [3].  We  may  write  17*  = e~tA,  where  —A  is 
the  infinitesimal  generator  of  17*.  Here,  A corresponds  to  the  fractional  differential 
operator  (— A)1//2.  Note  that  for  t > 0,  17*  maps  LP(R2)  into  D(A),  so  that  AUt  f is 
well-defined  for  arbitrary  f E Lp.  In  general,  this  is  not  the  case  for  non-holomorphic 
semigroups.  The  Gauss  singular  integral,  where  the  two-dimensional  Gaussian  kernel 
is  used  in  lieu  of  xb  in  (9),  defines  an  analogous  holomorphic  semigroup  IT*,  with 
A = —A.  Many  such  singular  integral  semigroups  5*  exist.  A very  rich  variety  can 
be  constructed  by  subordination  [7],  [40].  For  small  t > 0,  5*  behaves  as  an  approx- 
imate identity  on  Lp.  There  is  a large  literature  on  how  well  S*f  approximates  / 
as  t | 0.  See  [3],  [4],  [5],  [35],  [36],  [37],  and  the  references  therein.  As  t [ 0,  we 
have  ||  S*f  — / ||p=  o(l)  for  arbitrary  f E Lp , ||  Stf  — f ||p=  0(t)  if  and  only  if 

f E D(A),  and  ||  St f — f ||p=  o(t)  if  and  only  if  Stf  = f for  all  t > 0.  Thus, 
the  optimal  rate  is  always  0(t).  Of  particular  interest  in  this  paper  is  the  case  of 
non  optimal  approximation,  where  f </  D(A)  yet  retains  sufficient  smoothness  that 
||  Stf  — f ||p=  0(ta),  0 < q < 1.  as  t i 0.  While  complete  theories  exist  for  a wide 
class  of  singular  kernels,  the  simplest  such  theory  revolves  around  the  Poisson  semi- 
group Ut  in  (11).  We  have  from  [37,  Theorem  4], 

Theorem  1.  Let  17*,  t > 0.  be  the  Poisson  integral  operator  in  (11),  and  let 
0 < a < 1,  1 < p,  q < oo.  Then,  f E A (a,p,  q)  if  and  only  if 

rOC 

(12)  / (t~a  WlJtf  - f \\p)q  dt/t  < oo. 

Jo 

For  q = oo,  we  have  f E A (»,p,  oo)  if  and  only  if 

(13)  sup  t~a  ||  17*/  - / ||p  < oo. 

t>o 

Using  the  embedding  results  in  (7)  together  with  (13)  leads  to  the  following  corollary. 

Theorem  2 (Corollary).  Let  f E A (a,p,q),  with  2/(1  + a)  < p <2,  and  let 
(3  — 1 + a — p/2.  Then,  in  the  L2  norm 

(14)  . sup  t~ 3 ||  17*/  — / || 2 < oo. 

t>  o 


4.  Periodized  problems,  the  Poisson  summation  formula,  and  FFT  al- 
gorithms. The  above  results  can  be  used  to  fashion  a practical  image  analysis  tool. 
Theoretically,  given  any  image  f(x,y)  in  L1(R 2),  one  could  use  the  Fourier  transform 
(8)  to  form 

(15)  T {U*/}  = e-**7(£,  rf),  P = V^  + V2, 

for  sequences  of  positive  1- values  tending  to  zero.  Inverse  transformation  is  always 
possible  on  account  of  the  factor  e~tp , and  this  can  be  used  to  produce  an  infinite 
sequence  of  positive  numbers  pn  = {||  17*"/  — / ||  i / ||  / ||  1}  with  tn  [ 0.  If  every 
such  sequence  (tn,g.n),  ultimately  lies  below  the  curve  pit)  = C ta,  0 < t < t,  for 
suitably  chosen  constants  C > 0 and  0 < a < 1,  then  ||  17*/  — / ||i<  C ||  / ||i  ta , as 
t i 0,  and  f(x,y)  E A(q,  1,oo)  by  Theorem  1.  However,  this  does  not  lead  to  a 
practical  procedure. 
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On  the  other  hand.  Theorems  1 and  2 remain  valid  in  the  periodic  case  [36],  [37]. 
Here,  the  image  f(x,y)  and  the  kernel  i. jj(x,y,t)  in  (9)  are  now  periodized  [5],  [6].  Let 
Cl  denote  the  unit  square  — 1/2  < x,  y < 1/2  in  R2.  The  image  f(x,  y)  is  now  viewed 
as  originally  defined  on  Cl  from  which  it  is  extended  by  periodicity  to  all  of  R2.  Let 

(16)  f(Z,v)=  I f(x,y)e-2^x+™Uxdy. 

Jq 

Define  the  periodized  Poisson  kernel  ip*(x,y,t)  by 

OO 

(17)  ip*{x,y,t)=  Y x + k,y  + m,t ),  t > 0,  (x,y)eR2, 

k,m=  — oo 


and  let 

(18)  Utf=  / if* (x,y,t)f(x  — u,y  — v)dudv,  t > 0. 

J Q 

The  Poisson  summation  formula,  [1],  [5],  [6],  [22],  [37],  is  used  to  show  that  the 
periodized  Poisson  kernel  has  a complex  Fourier  series  with  Fourier  coefficients  again 
given  by  (10),  but  where  £,77  are  nowT  integers  running  from  —00  to  +00.  Moreover, 

OC 

(19)  Utf=  Y,  e-^/0e,r 7)e2«0*+OT),  t>  0,  p^^Y+V2- 

£,T]=-oo 

Again  the  factor  e~ip  assures  uniform  convergence  of  the  Fourier  series  in  (19).  Let 

N 

(20)  fN(x,y)  = Y e-tpna,v)e27ri(x(i+y7i\  t > 0,  p=  VC'  + V2- 

Z,V=~N 

Since  Lp( Cl)  C Zd(Q),  p > 1,  we  may  apply  this  approach  to  any  / e Lp,  and 
||  Ut f — /n  ||p  can  be  made  arbitrarily  small  by  choosing  N large  enough  in  (20). 
Next,  given  the  2 J x 2 J digitized  image  f(x,  y)  with  J > N,  the  discrete  Fourier  trans- 
form [2]  is  now  the  appropriate  numerical  tool  for  analyzing  this  periodized  problem. 
One  can  use  FFT  algorithms  to  form  the  Fourier  coefficients  /(£,  77),  — J < £,77  < J, 
and  then  apply  the  filter  ( e~tp  — 1 ) as  in  (15).  An  inverse  FFT  then  yields  an  accurate 
approximation  to  Ut f — f at  each  of  the  2J  x 2J  pixels,  for  each  small  t > 0.  We 
may  then  examine  the  discrete  Lp  relative  error  in  Poisson  approximation  as  t J.  0, 
and  locate  constants  C and  a such  that  ||  Utf  — f ||p<  C \\  f \\p  ta,  0 < t < t. 
In  summary,  we  have  constructed  an  accurate  numerical  procedure,  based  on  correct 
mathematical  analysis,  for  assessing  membership  in  any  A(a,p,  00)  space.  Equally 
important,  the  values  of  C and  a constitute  a-  priori  information  that  will  be  useful 
in  stabilizing  the  ill-posed  deblurring  problem. 

Remark  2.  Analogously  to  the  case  of  lossy  wavelet  compression  discussed  in  Remark 
1,  there  is  a finite  dimensionality  pitfall  in  the  above  singular  integral  methodology 
that  necessitates  the  exclusion  of  very  small  values  oft>  0.  Let  f°°(x,y ) be  the  orig- 
inal image  intensity  field  as  in  Remark  1,  and  assume  that  f°°(x,  y)  € A(0.5,p,  00),  so 
that  ||  £/f/°°  — /,oc  ||p=  O(Vt)  as  t \ 0,  by  Theorem  1.  Let  f(x,y)  be  the  2 J x 2 J 
digitized  image  corresponding  to  f°°(x,  y).  We  shall  show  that  at  very  small  values  of 
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t > 0,  the  behavior  of  ||  Utf  — f ||p  diverges  from  true  behavior  in  y ),  resulting 

in  a false  reading  for  a.  Let  St  = e~tA  be  any  contraction  semigroup  on  LP(R2).  As 
already  pointed  out,  if  / € D(A),  ||  S1  f — f ||p=  Oft)  as  t J,  0.  This  follows  from 

(21)  Stf-f=  [t-f(Suf)du=  [ SuAfdu, 

Jo  du  J o 

so  that  ||  St f — f ||p<  t ||  Af  ||p,  for  all  t > 0.  In  addition,  ||  S*f  — f ||p~  t ||  Af  ||p, 
for  all  sufficiently  small  t > 0,  because  ||  St  ||p~  1 for  all  small  t.  In  the  above  Poisson 
semigroup  t/*,  the  unbounded  operator  A is  defined  as  follows  in  Fourier  space 

(22)  r{Af}  = pfte,v),  p=ve  + v2- 

Since  the  digitized  2 J x 2 J image  f(x,  y ) is  a trigonometric  polynomial,  it  is  always 
€ D(A)  and  ||  Af  ||p  is  always  finite,  although  it  may  be  very  large.  Consequently, 
with  a possibly  large  positive  constant  K,  we  always  have  ||  Ui  f — f ||p<  Kt  for  all 
t > 0,  as  well  as  actual  linear  behavior  ||  \Jtf  — f ||p~  Kt  for  all  sufficiently  small 
t.  irrespective  of  the  behavior  of  ||  U1  f°°  — f°°  ||p  at  these  same  values  of  t.  This 
phenomenon  is  well-illustrated  in  Figure  1. 

5.  Application  to  real  images.  The  following  examples  illustrate  the  use  of 
the  Poisson  singular  integral  approach.  Our  first  example,  in  Figure  1,  is  the  512  x 
512  Mandrill  image  highlighted  in  [24]  as  an  example  of  an  image  ^ BV(R2).  The 
above  FFT  procedure  was  used  to  obtain  the  L 1 and  L2  relative  errors  in  Poisson 
approximation 

(23)  n(t)  =||  U*f-f  ||P/  ||  / ||P,  p — 1,2, 

at  300  values  of  t given  by  tn  = 0.5(0.95)n,  n — 1,300.  For  the  L 1 norm,  a plot 
of  pit)  versus  t on  a log-log  scale  produced  the  solid  curve  A in  Figure  1.  Least 
squares  fitting  was  used  to  find  the  two  distinct  majorizing  dashed  straight  lines  T 
and  E.  For  each  dashed  line,  the  y-axis  intercept  value  obtained  by  least  squares  was 
slightly  increased  so  as  to  make  each  line  lie  visibly  above  the  solid  curve  A\  however, 
the  slope  of  each  line  remains  the  same  as  that  obtained  from  least  squares.  The 
line  r,  defined  by  log  pit)  = 3.2  + 0.994  log  t,  accurately  captures  the  misleading 
linear  trend  of  p(t)  for  very  small  values  of  t,  while  being  grossly  inaccurate  at  larger 
values  of  t.  It  was  obtained  by  excluding  data  corresponding  to  log  t > —7  from 
the  least  squares  fit.  The  line  T implies  that  ||  U* f — f ||i<  24.53  ||  / ||i  t0-994 
for  all  t > 0.  As  emphasized  in  Remark  2,  this  correct  statement  primarily  reflects 
the  fact  that  the  512  x 512  Mandrill  image  lies  in  a finite  dimensional  space;  it  does 
not  describe  the  smoothness  properties  of  the  intensity  field  f°°(x,y)  that  gave  rise 
to  the  digitized  Mandrill  image.  The  majorizing  dashed  straight  line  E,  defined  by 
log p(t)  = — 0.75  + 0.306 log t,  accurately  reflects  behavior  for  —6  < logt  < — 1, 
while  being  grossly  inaccurate  at  very  small  values  of  t.  The  line  E was  obtained 
by  excluding  all  data  corresponding  to  logf  < —6  from  the  least  squares  fit.  Note 
that  this  still  leaves  over  100  data  points  remaining.  The  behavior  along  E,  where 
||  Lftf  — / ||i<  0.472  ||  / ||i  t0-306,  0 < t < 0.1,  is  taken  to  be  the  true  behavior  in  the 
Mandrill  image.  From  (13),  this  implies  that  the  Mandrill  image  € A(0.306, 1,  oo),  and 
hence,  is  not  of  bounded  variation.  The  behavior  in  the  L 2 norm  is  strikingly  similar, 
and  indicates  the  image  € A(0.271,  2,  oo).  Estimates  of  a in  any  other  discrete  Lp 
norm  can  be  obtained  similarly.  All  a estimates  shown  in  this  paper  were  obtained 
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Fig.  1.  512  x 512  Mandrill  image  was  identified  in  [24]  as  not  in  BV(R2).  This  is  confirmed 
in  above  graphical  use  of  Theorem  1,  using  FFT  techniques  discussed  in  Section  4-  Solid  curve  A is 
a plot  of  p(t)  =||  U*f  - f ||i  / ||  / ||i  versus  t,  on  a log-log  scale.  Majorizing  dashed  straight  line 
F,  defined  by  log  p(t)  = 3.2  + 0.994 log  t,  accurately  captures  linear  behavior  of  p(t)  for  very  small 
values  of  t,  but  is  grossly  inaccurate  at  larger  values  of  t.  Linear  behavior  at  very  small  values  of 
t is  misleading,  and  falsely  implies  that  true  image  is  of  bounded  variation.  (See  Remarks  1 and 
2).  Majorizing  dashed  straight  line  E,  defined  by  log  p(t)  = —0.75  + 0.306  log  t,  accurately  reflects 
behavior  for  —6  < log  t < —1,  while  being  grossly  inaccurate  at  very  small  values  oft.  Behavior 
along  E.  where  ||  Ut  f — f ||i<  0.472  ||  / ||i  t0306,  0 < t < 0.1,  is  taken  to  be  true  behavior  in 
Mandrill  image.  From  (13),  this  implies  image  € A(0.306, 1,  oo). 


using  the  above  procedure  of  constructing  the  line  E in  log-log  plots  of  ft(t),  after 
excluding  all  data  corresponding  to  logf  < —6.  As  in  [18,  §5B],  we  have  occasionally 
found  contradictory  examples  where  the  value  of  a in  the  L2  norm  was  greater  than 
that  in  the  Ll  norm.  When  that  happened,  a new  E line  was  constructed  for  the  L2 
trace,  based  on  excluding  data  corresponding  to  log  t < — 5.  It  is  recommended  that 
data  for  very  small  values  of  t always  be  included  in  log- log  plots  of  so  as  to 
enable  clear  identification  of  the  spurious  linear  trend  prior  to  rejecting  that  part  of 
the  data. 

Our  second  example,  in  Figure  2(A),  is  a 1024  x 1024  Whirlpool  Nebula  image, 
taken  at  the  National  Optical  Astronomy  Observatory,  (NOAO/AURA/NSF),  by  T. 
Rector  and  hi.  Ramirez.  As  in  the  case  of  Figure  1,  Poisson  integral  approximation 
in  Ll  was  used  to  obtain  the  solid  curve  A,  and  the  line  E^  was  constructed  using 
least  squares.  This  procedure  was  repreated  for  the  L2  norm.  The  results  indicate 
that  Figure  2(A)  satisfies  ||  Ut f — f ||i<  0.6  ||  f ||i  t °'530,  0 < t < 0.1,  and  that 
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Fig.  2.  Whirpool  Nebula  M51.  Original  and  sharpened  images  have  noticeably  different  L 1 
Poisson  traces  p(t)  = ||  U* f — / ||i  / ||  / ||i,  reflecting  substantially  different  values  of  Lipschitz 

exponent  a.  (A)  Original  1024  x 1024  image  taken  by  T.  Rector  and  M.  Ramirez,  National  Optical 
Astronomy  Observatory,  (NOAO/AURA/NSF).  Poisson  approximation  produces  solid  trace  A,  ma- 
jorized by  dashed  straight  line  E^  defined  by  \ogp(t)  = — 0.5  + 0.530  log  t.  This  implies  that  image 
(A)  € A(0.530, 1,  oo).  (B)  Blind  deconvolution  of  (A)  using  APEX  method  [9],  brings  out  signifi- 
cant fine  scale  detail,  and  results  in  solid  trace  B,  majorized  by  dashed  straight  line  E b defined  by 
log p(t)  — —0.2  + 0.239  log  t.  This  indicates  that  deblurred  image  (B)  € A(0.239,  l,oo).  Image  (B) 
strongly  resembles  [3f,  plate  26]  taken  by  Milton  Humason  using  200  inch  Mt  Palomar  telescope. 
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Fig.  3.  A significant  class  of  high-resolution  8-bit  images  have  Lipschitz  exponents  a in  the 
range  0.2  < a < 0.6,  in  either  L 1 or  L~ , and  are  not  of  bounded  variation. 


TABLE  1.  Values  of  (C,  a)  in  ||  Ut  f — f ||p<  C ||  / ||p  ta,  0 < t < 0.1,  p = 1,  2, 
for  each  image  f(x,y ) in  Figure  3,  when  U*  is  Poisson  operator  in  (19). 


Image 

Size 

(C, 

a)  e 

A(a, 

l,oo) 

(C, 

ct)  G 

A(a, 

2,oo) 

Marilyn  Monroe 

5122 

c = 

0.77, 

a = 

0.565 

c = 

0.68, 

a = 

0.474 

Sagittal  brain  MRI 

5122 

c = 

1.28, 

Q = 

0.590 

c = 

1.02, 

a = 

0.520 

Washington  DC  Landsat 

5122 

c = 

0.45, 

a = 

0.341 

c = 

0.55, 

a = 

0.340 

Mariner  5 spacecraft 

5122 

c = 

0.90, 

a = 

0.448 

c = 

0.99, 

a — 

0.417 

USS  Eisenhower 

5122 

c = 

0.47, 

a = 

0.420 

c = 

0.50, 

a — 

0.362 

English  Village 

5122 

c = 

0.49, 

a = 

0.472 

c = 

0.55, 

Q = 

0.439 

Nanoscale  micrograph 

10242 

c = 

0.45, 

a = 

0.415 

c = 

0.55, 

a = 

0.415 

Spiral  galaxy  MSI 

10242 

c = 

0.68, 

a = 

0.365 

c = 

0.78, 

a — 

0.327 

Cluster  of  galaxies 

10242 

c = 

0.65, 

a = 

0.222 

c = 

0.97, 

a = 

0.216 

SINGULAR  INTEGRALS  AND  RECOVERY  OF  IMAGE  TEXTURE 


11 


Figure  2(A)  € A(0.530, 1,  oo)  (T  A(0.462,  2,  oo).  Interestingly,  if  we  sharpen  Figure 
2(A)  using  the  APEX  method  [9],  we  obtain  the  image  in  Figure  2(B).  This  enhanced 
image  displays  significant  fine  scale  detail  not  readily  visible  in  the  original  image, 
and  strongly  resembles  a Whirlpool  galaxy  plate  taken  by  Milton  Humason  in  1950 
using  the  200  inch  Mt  Palomar  telescope.  See  [34,  plate  26).  Here,  L 1 Poisson  analysis 
produced  the  solid  curve  B and  the  majorizing  line  S b-  We  find  that  Figure  2(B) 
€ A(0.239, 1,  oo)  n A(0.230,  2,  oo),  and  thus  has  substantially  lower  values  of  a than 
does  Figure  2(A).  This  result  is  highly  plausible.  Presumably,  any  low-pass  blurring 
process  that  may  have  affected  Figure  2(A)  would  have  attenuated  fine  scale  features, 
and  thereby  artificially  increased  the  values  of  a.  The  result  also  indicates  that  APEX 
processing  of  image  (A)  produced  relatively  more  sharpening  in  the  Ll  norm  than  in 
the  L 2 norm. 

The  nine  images  in  Figure  3 and  Table  1 form  an  interesting  collection  that 
includes  natural  as  well  as  man  made  objects,  exhibiting  a wide  range  of  sizes.  The 
last  row  contains  an  electron  microscop}'  nanoscale  structure,  a galactic  scale  object, 
and  a cosmological  scale  structure.  Along  with  the  three  images  in  Figures  1 and 
2,  this  paper  has  applied  the  Poisson  integral  method  to  12  high  resolution  images, 
and  we  have  found  that  in  either  A(a,  1,  oo)  or  A(a,  2,  oo),  the  values  of  a lie  in  the 
range  0.2  < a < 0.6.  This  range  of  values  is  compatible  with  that  found  in  [12], 
[18],  [19],  using  an  entirely  different  method.  Moreover,  while  A(a,2,oo)  are  smaller 
spaces  than  are  A(ct,  l,oo),  they  are  evidently  wide  enough  to  contain  each  of  these 
12  images,  albeit  with  smaller  values  of  a.  Notice  also  that  the  values  of  the  constant 
C in  Table  1 are  confined  to  a very  narrow  range  in  both  Ll  and  L1 . 

The  PSI  deblurring  method  to  be  described  in  Section  7 below  requires  prior 
knowledge  of  the  values  of  C and  a in  the  desired  unknown  deblurred  image.  In 
general,  these  values  will  not  be  known.  However,  as  shown  in  Figure  2,  it  is  reasonable 
to  assume  that  a in  the  deblurred  image  will  be  lower  than  in  the  given  blurred 
image,  provided  that  image  is  relatively  noise  free.  Moreover,  given  a fast  algorithm, 
because  of  the  narrow  range  of  values  involved  in  both  C and  a,  it  is  feasible  to 
simultaneously  compute  and  display  multiple  reconstructions,  based  on  numerous 
hypothetical  choices  of  C'  and  a.  Efficient  exploration  in  parameter  space  is  often  the 
key  to  the  successful  solution  of  inverse  problems,  when  such  problems  can  be  solved. 

6.  Image  deblurring  in  L2(R2).  We  now  consider  the  image  deconvolution 
problem  Pf  = g with  a known  shift-invariant  point  spread  function  (psf)  p(x,y ), 

(24)  Pf  = p(x,y)  <g>  f(x,y)  = g(x,y),  g(x,y)  = ge{x,y)  + n{x,y). 

Here,  <g>  denotes  convolution,  g(x,  y)  is  the  given  recorded  noisy  blurred  image,  ge(x,  y) 
is  the  hypothetical  exact  blurred  image  that  would  have  been  recorded  in  the  absence 
of  any  noise,  and  n(x , y),  presumed  small,  represents  the  cumulative  effects  of  all  noise 
proceses  and  other  errors  affecting  final  acquisition  of  the  digitized  array  g(x,  y).  The 
noise  may  be  multiplicative.  Neither  ge(x,y)  nor  n(x,y)  are  known,  only  their  sum 
g(x,  y).  Denoting  the  unknown  exact  sharp  image  by  fe(x,  y),  we  have 

(25)  Pfe  = p(x,  y)  <g>  fe(x,  y)  = ge(x,  y). 

Given  only  (24),  we  seek  a solution  f(x,  y)  in  (24)  such  that  Pf  s;  g , and  such  that 
II  / — fe  ||  2 is  small.  To  achieve  this  goal,  some  a-priori  information  about  fe  and 
n is  always  necessary.  Most  real  images  fe(x,y)  contain  fine-scale  features,  sharp 
edges,  and  other  kinds  of  non-differentiable  singularities.  Deblurring  techniques  that 
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impose  stabilizing  constraints  in  the  form  of  prescribed  bounds  on  partial  derivatives 
of  f(x,y ) in  (24),  are  generally  inapplicable,  although  they  are  often  used.  Penalties 
for  such  use  include  smoothing  out  of  sharp  features,  and  possible  loss  of  vital  diag- 
nostic information.  Indeed,  the  desire  to  accurately  reconstruct  edges  and  other  sharp 
singularities  was  the  principal  reason  for  developing  total  variation  methods.  In  fact, 
several  deblurring  methods  actually  exist  that  do  not  require  bounds  on  derivatives 
[8]. 

We  shall  limit  discussion  of  (24)  to  the  case  of  uniform  defocus  blur,  where  the 
psf  is  proportional  to  the  characteristic  function  of  a disc  of  radius  R.  This  is  the 
so-called  ‘pillbox’  model  [17],  [25],  [16],  [31].  If  R > 0 is  the  radius  of  the  ‘circle  of 
confusion’,  the  psf  for  defocus  blur  is  given  by 

( (7ri?2)-1,  x2  + y2<R2, 

(26)  P0r,y)=< 

i 0,  x2  + y2>R2. 

This  has  a Fourier  transform  given  by  the  ‘sombrero  function’  [23,  p.  72] 

(27)  p(Z,  77)  = 2Ji  (Rp)/(Rp),  p = v^2  + V2, 


where  Ji(x)  is  the  Bessel  function  of  the  first  kind  of  order  1.  In  our  numerical 
experiments  below  on  2N  x 2 N images,  the  expression  (27)  is  used  to  blur  images  by 
Fourier  domain  multiplication  with  a preselected  R > 0,  and  £,  77,  are  integers  with 
-N  < £,77  < N.  Rather  than  interpret  R as  a radius,  we  simply  observe  that  the 
severity  of  such  a blur  is  determined  by  the  number  of  zeroes1  in  \p(p)\  on  0 < p < N . 

6.1.  True  Wiener  filtering  and  the  Tikhonov- Miller  method.  Wiener 
filtering  [32,  p.  356],  is  an  important  example  of  a method  that  does  not  impose 
differentiability  constraints.  It  assumes  instead  that  the  power  spectra  |n(£,  77)  | and 
l/e (C'7?) I °f  each  of  n(x,  y)  and  fe(x,y)  are  known.  When  this  is  the  case,  Wiener 
filtering  produces  a solution  fw(x,  y)  in  (24)  defined  as  follows  in  Fourier  space 


(28) 


r&v) 


p(^v)g(^v) 

\p(t,V)\2  + \n(t,v)\2/\fe(t,v)\2 


where  3 denotes  the  complex  conjugate  of  2.  Under  some  additional  conditions,  it 
can  be  shown  that  fw(x,  y ) is  an  approximate  solution  of  Pf  — g that  minimizes 
the  error  ||  / — fe  H2  over  all  / € L2.  In  practice,  the  power  spectra  |n(£,  77)!  and 
|/e(£,  77)1  are  very  seldom  known  in  advance,  and  true  Wiener  filtering  is  almost  never 
realizable.  However,  the  solution  (28)  is  of  considerable  theoretical  interest  because 
of  its  optimality  property.  Note  that  numerous  ad-hoc  versions  of  (28)  exist,  in  which 
more  readily  available  quantities  are  substituted  in  place  of  the  required  but  unavail- 
able true  power  spectra.  Such  versions  are  sometimes  called  Wiener  filtering  by  an 
abuse  of  terminology.  However,  these  substitute  versions  do  not  satisfy  the  Wiener 
optimality  criterion,  nor  do  they  elicit  the  same  degree  of  theoretical  interest. 

One  of  the  best-known  rigorously  analyzable  and  feasible  versions  of  Wiener  filter- 
ing is  the  Tikhonov-Miller  method  [28],  now  considered  canonical  in  image  deblurring 
[25] . Significantly,  this  method  makes  no  a-priori  assumptions  regarding  the  statistical 
character  of  the  data  noise.  For  non-differentiable  images,  Tikhonov-Miller  restora- 
tion requires  the  following  a-priori  information:  an  upper  bound  e > 0 for  the  L2 


1The  first  five  positive  zeroes  of  J\(x)  are  3.83171,  7.01559,  10.17347,  13.32369,  and  16.47063. 


SINGULAR  INTEGRALS  AND  RECOVERY  OF  IMAGE  TEXTURE 


13 


norm  of  the  noise  n(x,  y)  in  the  blurred  image  g(x,  y),  and  an  upper  bound  M for  the 
L 2 norm  of  the  unblurred  image  fe 

(29)  ||  n ||2=||  pfe  - 9 lb<  G ||  fe  ||2<  M,  e/M  < 1. 

It  is  assumed  that  e and  M are  compatible  with  the  existence  of  an  fe(x,y)  € L2 
satisfying  (29).  Tikhonov-Miller  restoration  is  defined  as  the  unique  function  fT(x,  y) 
such  that 

(30)  fT(x,  y)  = Arg  min  { ||  Pf  - g \\l  +(e/M)2  ||  / ||1}  . 

j (zL-(R- ) 

As  will  be  seen  from  Theorem  3,  where  the  Tikhonov-Miller  method  corresponds  to 
the  special  case  Tj-  = 0,  this  minimum  problem  has  a unique  solution  satisfying 

(31 ) QrfT  = P*g,  Qr  = P*P  + (e/M)2/. 

Moreover,  there  holds  the  following  best-possible  error  bound  for  Tikhonov-Miller 
reconstruction 

(32)  II  fT  - fe  ||2<  ex/2  II  Qf1/2  ||2, 


where 

(33)  \\Qf1/2  \\2=swp{\p(£,ri)\2 + {e/M)2}  1/_. 

C»7 


Given  the  psf  p(x,  y),  together  with  the  a-priori  information  e,  M,  one  can  always  find 
the  maximum  value  in  the  2 N x 2 N array  on  the  right  of  (33).  As  in  (28)  we  may 
implement  (31)  in  Fourier  space.  We  have 


(34) 


r&v) 


P( ^ V)9(Z,  9) 

|p(e,fy|2  + (e/A/)2- 


Moreo\’er.  from  (29)  and  Parseval's  relation 

(35)  [ |fi(£)  g)\2dfdg  < e2 , [ \fe(£,  g)\2dfdg  < M2. 

Jr - Jr 2 

Therefore,  the  Tikhonov-Miller  method  can  be  viewed  as  an  approximate  version  of 
true  Wiener  filtering  where  the  unavailable  pointwise  values  of  the  spectra  in  (28)  are 
replaced  by  more  readily  available  integrals  of  these  quantities.  However,  it  may  be 
anticipated  that  since  true  Wiener  filtering  requires  prior  knowledge  in  the  form  of 
8N2  numbers  for  a 2N  x 2N  image,  whereas  the  Tikhonov-Miller  method  requires 
only  2.  less  accurate  results  must  generally  be  expected  from  the  latter  method. 

7.  The  Poisson  Singular  Integral  (PSI)  method  for  images  e A(q,  2,  00). 
The  preceding  discussion  was  necessary  to  set  the  stage  for  the  PSI  method.  Here,  in 
addition  to  the  a-priori  constraints  (29),  the  behavior  of  ||  U*  fe  — fe  ||2  on  0 < t < t 
is  assumed  known,  as  in  the  case  of  Table  1 . The  constants  Cj  and  a are  now  used 
to  place  a further  constraint  on  fe(x,  y).  For  any  / € L2(R2),  we  have  on  Fourier 
transforming  / — Uf  f and  using  (10), 


(36) 


JF  {/  - t/fy}  = (1  - e~tp)  f(t  V),  P = Ve  + P2- 
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I R- 


Therefore,  from  Parseval’s  theorem, 

(37)  [ ||  Usf  - f \\l  ds  = [ ds 

Jo  Jo 

For  fixed  t > 0,  define  £(£,  77,  t)  > 0 by 

(38)  z(^,rj,t)  = \[  (l  — e~sp)~  ds 


(1 


-sp 


Y \f(Z,v)\2d£dri. 


t + 


4e~tp  — e~2tp  — 3 ) 
2 P j 


1/2 


It  follows  directly  from  the  integral  definition  in  (38)  that  for  any  fixed  t > 0,  z(p,  t) 
is  a strictly  increasing  function  of  p,  and  that  i(0,0,  t)  — 0.  For  fixed  t > 0,  define 
the  linear  operator  Z(t)  in  L2(R2)  by 

(39)  Z(t)f=  [ z(^V,t)f(^V)e2^x+^d^dV. 

Jr- 

Then,  from  (37), 


(40) 


Us f — f || 2 ds  =||  Z(t)f  ||i. 


For  any  fe  & A(q,  2,  00).  0 < a < 1,  we  have  ||  Usfe  — fe  ||2<  Cj  ||  fe  H2  sQ,  0 < s < t, 
where  Cj  is  a positive  constant  depending  on  t,  fe  and  a.  Therefore,  with  ||  fe  ||2<  M 


>9  71  .r9Tl  + 2a 


(41) 

\\Z(t)fe\\l<  t 

II  \ jje  112-  1+2ft 

Define 

(42) 

( \ V2 

r-=  / 1 + 2q  1 
\ c=Y+2a  J 

Then, 

(43)  ||  P/e  - g ||2<  e, 

(e/M)  ||  /e  ||2<  e,  l 

FLx  t > 0,  and  consider  the  minimization  problem 

(44)  P'(x,y)  = Arg  min  o {\\Pf-g  |||  + (e/7\/)2  (||  f \\%  +r|  ||  Z(t)f  |||) } - 

f GL-(R-) 

As  will  be  seen  in  Theorem  3 below,  this  minimum  problem  has  a unique  solution 
satisfving 


(45) 


Qpf  = p*3,  Q*  = P*P  + (e/A/)2  {/  + r=z(iyz(i)} . 


The  function  f^(x,y)  in  (44)  is  defined  to  be  the  PSI  deblurred  image.  Moreover, 
there  holds  the  following  error  bound  for  PSI  deblurring 


(46) 

where 


f*  - fe  l|2<  ev/3  j|  <j:1/2 


t\,2\  1-1/2 


Qw  ' l|2=  sup  {|p(C,7?)|2  + (e/fi/)2  (1  + r=|c(^,7?,t)|2)} 


(47) 
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Given  the  psf  p(x,y),  together  with  the  a-priori  information  e,  M,  and  Tj,  one  can 
always  find  the  maximum  value  in  the  2N  x 2N  array  on  the  right  of  (47).  Again,  as 
in  (28)  and  (34),  can  be  found  explicitly  in  Fourier  space.  We  have 


(48) 


Hit  n\  = flMj) 

1 m,v)\2  + (e/M)2{!  + rf|£(e,r?J)|2}’ 


Equations  (45-48)  should  be  compared  with  equations  (31-34).  Tikhonov-Miller  de- 
blurring  can  then  be  seen  as  an  extreme  case  of  PSI  deblurring,  the  case  where  fe(x,  y) 
is  presumed  no  smoother  than  the  most  general  L1  function,  so  that  ||  Utfe  — fe  || o= 
o(l)  as  t { 0.  This  corresponds  to  Cj  = oo  in  (41),  and  hence,  Tj  = 0 in  (48). 


Theorem  3.  Fix  t > 0 and  let  the  exact  image  fe(x,y)  satisfy  the  a-priori 
constraints  (43).  Let  f^'(x,y)  minimize  (44)>  and  let  Qq,  be  the  positive  self-adjoint 
operator  on  L2(R2)  given  by 

(49)  Qj,  = P*P+  (e/M)2  {I  + Tf  Z(?)*Z(i)}  . 

Then  fP  is  the  unique  solution  of  Qq,f^  = P*g,  and  satisfies 

II  Pf  - 9 III  +(e/M)2  { ||  ft  III  +T|  ||  zq£)f+  III}  < 3e2, 

(50) 

II  Pif  ~ fe)  III  +(e/M)2  { ||  /*  - fe  II2  +Tf  ||  Z(t)(f*  - fe)  ||2}  < 3e2. 

This  implies  the  L 2 error  bound 

(51)  ||  - fe  ||2<  e v/3  ||  Qf1’2  ||2, 

where 

(52)  ||  Qfl/2  ||  2=  sup  { [p(£,  17)  |2  + (e/M)2  (l  + r5|£(^,  77,  i)|2)  }_1/2  . 

€,v 

Proof.  Let  H denote  the  Hilbert  space  direct  sum  L2(R 2)  0 L2(R2)  0 L2(R2) 
with  elements  [u,v,w\,  scalar  product  ([tti,  v\,  tui],  ]v,2,  V2,  wff)  =<  u\,uo  > + < 
v\,V2  > + < u>i , W2  >,  and  norm  |||  |||.  Let  P : L2(R2)  1— > TL  be  defined  by 

Pf  = [Pf,  iwf,  uTjZ(t)f],  where  uj  — (e/M),  and  let  g = [p,  0, 0].  We  seek  to 
minimize  |||P/  — g|||  over  all  / G L2(R2).  The  normal  equation  P* Pf^  = P*g  gives 
Qii’f^  — P*9  with  Q as  in  (49).  By  hypothesis  \\\Pfe  — g|||2  < 3e2.  The  minimizing 
element  is  such  that  Pf^’  is  the  orthogonal  projection  in  Tt  of  g on  the  range  of 
P.  By  the  Pythagorean  theorem 

(53)  |||P/*  - 5|||2  + III P(fe  - /*)|||2  = 1 1 1 P fe  ~ 5|||2  < 3e2. 

This  proves  (50).  We  now  establish  (51).  From  (49),  (50), 

(54)  II  (?72(/e  - f*)  ||1=<  <?*(/,  - f+),  (u  - f*)  >=  |||P(/«  - /Pill2  < 3e2. 


Hence, 
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Fig.  4.  Instructive  deblurring  experiment  with  exact  a-priori  information  highlights  significant 
differences  in  behavior  in  above  three  FFT-based  methods.  (A)  Defocused  Marilyn  Monroe  image 
with  R = 0.06  and  3%  multiplicative  noise.  (B)  Tikhonov- Miller  method  with  exact  parameters  e and 
M.  brings  out  significant  noise.  (C)  PSI  method  with  exact  parameters  e,  M,  a = 0.474,  Cj  = 0.68. 
(D)  True  Wiener  filtering  with  exact  power  spectra  |n(£,  rj)|,  |/e(£,i?)|-  Realizable  PSI  deblurring 
closely  matches  unrealizable  true  Wiener  filtering. 

TABLE  2. 

Behavior  in  defocused  Marilyn  Monroe  image  in  Figure  4 ■ 


Deblurring  Method 

L 1 relative  error 

L 1 relative  error 

Tikhonov-Miller  (B) 

29.82% 

34.17% 

Poisson  Singular  Integral  (C) 

6.81% 

9.03% 

True  Wiener  Filtering  (D) 

6.03% 

7.88% 
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ERE  OR  BOUND  FUNCTIONS  ( EQ . 56)  IN  FIGURE  4 EXPERIMENT. 


DISCRETE  FREQUENCY  ( INTEGER ) 

POISSON  SINGULAR  {SOLID)  TIKHON  0 V— MILLER  ( DASHED ) 


Fig.  5.  Plot  of  error  bound  functions  0T(f),  (dashed  curve),  and  (solid  curve),  as  defined 

in  Eq.  (56),  for  the  deblurring  experiment  in  Figure  4 ■ Maximum  value  in  6T  is  more  than  six  times 
larger  than  in  6 ^ . Qualitative  difference  in  behavior  in  these  two  curves  implies  significant  difference 
in  Fourier  domain  regularization  in  the  PS  I and  Tikhonov- Miller  methods.  Difference  in  maximum 
values  explains  large  difference  in  L 2 relative  errors  in  Figures  4(B)  and  4(C). 


8.  A preliminary  deblurring  experiment.  In  the  following  controlled  exper- 
iment, knowledge  of  the  exact  solution  fe(x,y)  is  used  to  derive  exact  values  for  all 
parameters  that  constitute  a-priori  information  in  each  of  the  above  three  methods. 
Such  exact  knowledge  is  not  available  in  practice.  The  experiment  is  primarily  of 
theoretical  interest.  It  is  designed  to  illustrate  major  differences  in  behavior,  and  to 
properly  locate  the  PSI  method  in  relation  to  Wiener  filtering  and  the  Tikhonov-Miller 
method. 

The  8-bit  512  x 512  Marilyn  Monroe  image  fe{x,y)  in  Figure  3 was  synthetically 
defocused  by  Fourier  domain  multiplication  with  the  expression  in  (27)  using  R = 0.06. 
This  produced  the  exact  blurred  image  ge(x , y).  Multiplicative  noise  n(x,  y)  was  then 
added  to  ge  as  follows.  Each  pixel  value  ge(x,y ) was  perturbed  by  adding  to  it 
the  quantity  n(x,y)  = 0.03cr(x,y)ge(x,y),  where  a(x,y)  is  an  array  of  uniformly 
distributed  random  numbers  in  the  range  [—1,1].  We  term  this  process  ‘adding  3% 
noise’.  With  varying  percentages,  we  shall  use  the  same  process  in  all  our  experiments. 
Note  that  no  noise  is  thereby  added  at  points  where  ge(x,y)  = 0.  The  resulting 
g(x,y ) = ge(x,y ) + n(x,y ) is  shown  in  Figure  4(A).  We  find  ||  n || 2=  e = 2.247,  and 
||  fe  ||  1=  107.59,  ||  /e  || 0=  M = 131.13.  Therefore  e/M  — 0.01713.  From  Table 
1,  we  have  ||  Ulfe  — fe  || 2 < 0.68  ||  / H2  t0'4'4,  0 < t < 0.1.  With  t = 0.1,  (42) 
gives  Tj  = 19.33.  Next,  using  FFT  algorithms,  we  obtain  the  exact  power  spectra 
|n(£,p)|,  |/e(£,p)|.  We  are  now  ready  to  compare  these  three  FFT-based  procedures 
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under  optimal  conditions  for  each  method. 

The  Tikhonov-Miller  reconstruction  is  shown  in  Figure  4(B).  Significantly,  this 
reconstruction  is  quite  noisy,  despite  the  use  of  exact  values  for  e and  M.  While  the 
regularizing  information  in  (29)  prevents  explosive  noise  amplification,  it  is  obviously 
insufficient  to  prevent  serious  noise  contamination.  This  is  generally  the  case  in  the 
Tikhonov-Miller  method.  The  Poisson  Singular  Integral  restoration  is  shown  in  Figure 
4(C).  Here,  the  additional  information  that  fe  E A(a,  2,  oo),  together  with  the  values 
of  the  constants  Cj  and  a,  were  evidently  decisive  in  eliminating  noise.  The  Wiener 
filtered  solution,  shown  in  Figure  4(D),  appears  only  slightly  better  than  the  PSI 
solution.  However,  the  very  major  difference  between  true  Wiener  filtering  and  the 
approximate  version  known  as  the  Tikhonov-Miller  method,  is  another  significant 
result  brought  out  by  this  deblurring  experiment. 

It  is  instructive  to  study  the  Ll  and  L2  relative  error  pattern  shown  in  Table  2. 
It  is  widely  assumed  in  practice  that  the  L2  minimum  error  property  of  true  Wiener 
filtering  remains  more  or  less  true  for  the  more  feasible,  approximate  versions  of  such 
filtering.  This  is  not  the  case.  The  Tikhonov-Miller  relative  errors  are  more  than  four 
times  larger  than  the  true  Wiener  errors.  On  the  other  hand,  relative  errors  in  the 
PSI  method  are  only  slightly  larger  than  those  for  true  Wiener  filtering.  Put  another 
way,  the  PSI  method  appears  to  be  a feasible  procedure  that  can  very  substantially 
improve  upon  the  Tikhonov-Miller  method. 

Insight  into  how  this  improvement  comes  about  can  be  gained  by  an  analysis  of 
the  respective  error  bounds  for  each  method.  Notice  that  each  of  the  denominators 
on  the  right  hand  sides  of  (34)  and  (48)  are  radially  symmetric  functions  of  (£,  77), 
while  this  is  not  the  case  in  (28).  These  denominators  play  a dual  role.  They  define 
the  actual  regularization  procedures  in  (34)  and  (48),  and  they  define  the  resulting 
error  bounds  in  (33)  and  (47).  Because  of  the  radial  symmetry,  a one-dimensional 
picture  tells  the  whole  story.  Define  the  respective  Tikhonov-Miller  and  PSI  error 
bound  functions  9T(£),  9^ (£)  as  follows 


9r(Z)  = {mO)\2  + (e/M)2}  1/2  . 


(56) 


In  Figure  5,  we  plot  9r(f)  and  9U, ( f ) as  determined  by  the  actual  parameter  val- 
ues that  entered  the  deblurring  experiment  in  Figure  4.  The  significant  differences 
in  these  two  curves  translate  into  fundamental  differences  in  the  Fourier  space  reg- 
ularization that  defines  the  corresponding  procedures.  From  (27),  we  see  that  #r(£) 
has  a maximum  of  M/e  = 58.36,  at  every  point  £ > 0 where  Ji(0.06  f)  = 0.  There 
are  4 such  points  on  0 < £ < 256.  The  curve  #p(£)  also  develops  maxima  at  these 
same  points,  but  these  maxima  are  about  six  times  smaller  than  those  in  #r(£),  owing 
to  the  additional  term  involving  z(£,  0,  t).  Since  the  error  estimate  in  each  method 
is  proportional  to  the  maximum  along  the  corresponding  curve,  it  is  natural  to  find 
substantially  smaller  errors  in  Figure  4(C)  than  in  Figure  4(B). 

9.  Comparing  total  variation  deblurring  with  PSI  deblurring.  The  use  of 

initial  value  PDE  methods  in  image  processing  and  computer  vision  has  mushroomed 
into  an  important  new  branch  of  applied  mathematics.  The  basic  idea  originates  in 
gradient  descent  methods  for  minimizing  appropriate  energy  functionals.  Instructive 
surveys  of  this  general  set  of  ideas  may  be  found  in  [11]  and  [39], 
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Fig.  6.  Comparison  of  total  variation  and  PSI  deblurring  on  mildly  blurred  image.  Zooming  on 
selected  parts  of  the  image  enables  meaningful  comparisons  between  the  two  methods.  (A)  Original 
sharp  USS  Eisenhower  image.  (B)  Mildly  defocused  image  with  R = 0.03  and  0.1%  multiplicative 
noise.  (C)  Total  variation  deblurring  of  image  (B)  by  applying  finite  difference  scheme  in  [26,  §5/ 
to  Eq.  (59),  with  0 = 0.0001,  A = 300,  At  = O.l(Ax)2,  T = lOOAl.  (D)  PSI  deblurring  of  image 
(B)  using  exact  a-prion  parameters,  e,  M,  a = 0.362,  Cj  = 0.50.  (E)  Zooming  in  TV  deblurred 
image  reveals  significant  loss  of  structural  detail.  (F)  Zooming  on  same  region  in  PSI  image. 
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Fig.  7.  Comparison  of  total  variation  and  PSI  methods  on  moderately  blurred  image.  Zooming 
now  reveals  unacceptable  loss  of  content  in  TV  deblurring.  (A)  Original  sharp  English  Village  image 
(E)  Moderately  defocused  image  with  R = 0.06  and  0.1%  multiplicative  noise.  (C)  Total  variation 
deblurring  of  image  (B)  using  scheme  m [26,  §5/  with  0 = 0.0001,  A = 500,  At  = 0.1(Ax2),  T = 
lOOAt.  (D)  PSI  method  with  exact  a-priori  parameters,  e,  M,  a = 0.439,  Cj  = 0.55  (E)  Zooming 
in  image  (C)  reveals  loss  of  windows  and  roof  shingles.  (F)  Zooming  on  same  region  in  PSI  image. 
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The  total  variation  approach  introduced  in  [33]  is  one  of  the  most  popular  PDE 
methods,  and  it  is  primarily  designed  to  recover  edges  in  the  original  image.  Given 
the  deconvolution  problem  Pf  = g as  in  (24),  TV  deblurring  presupposes  the  exact 
sharp  image  fe(x,y ) G BV(R2),  and  it  produces  an  Image  ftv(x,y)  defined  by 


(57)  ftv(x,  y)  = Arg  min  { (A/2)  ||  Pf  - g \\%  + f I Vf\dxdyX  . 

feBV(Ri)  ( JR 2 J 

This  means  that  ftv(x,y)  is  the  solution  of 

(58)  P'Pf*  - A-'V.  (j!(Z|)  =P'9 


Here,  A > 0 is  a regularization  parameter  that  can  be  tuned.  Provided  the  noise 
level  in  g(x,  y)  is  small,  larger  values  of  A produce  sharper  images.  Too  large  a value 
of  A leads  to  computational  instability.  Unlike  the  cases  in  (31)  and  (45),  (58)  is  a 
nonlinear  deconvolution  problem  that  cannot  be  solved  explicitly  in  Fourier  space. 
In  fact,  considerable  effort  is  generally  required  to  obtain  ftv  for  large  size  imagery. 
In  p\:re  denoising  applications,  where  P — /,  this  effort  is  usually  warranted  by  the 
quality  of  the  resulting  restoration.  Recent!}',  a new  time  dependent  evolutionary 
approach  to  (58)  has  been  developed  [26],  whereby  ftv(x , y)  is  obtained  as  the  steady 
state  solution  to  the  following  nonlinear  anisotropic  diffusion  problem 


(59) 


' ut  = — A|Vw|  P*(Pu  -g)  + |Vw|  V.  (Vu/{  VlVu]2  + /?})  , 

< 

k u(x,y,  0)  = g(x,y), 


where  the  given  blurred  image  g(x,  y)  is  used  as  the  initial  value.  In  addition,  u(x,  y,  t ) 
satisfies  homogeneous  Neumann  conditions  at  the  boundary  of  the  unit  square  Q.  In 
(59),  (3  > 0 is  a small  constant  designed  to  prevent  division  by  zero.  In  [26,  §5],  an 
efficient  new  explicit  finite  difference  scheme  for  (59)  is  proposed.  This  scheme  has 
improved  stabilty  and  edge-enhancing  properties,  and  converges  rapidly  to  the  desired 
steady  state  solution.  Accordingly,  we  shall  use  that  method  in  our  total  variation 
deblurring  experiments. 

This  paper  has  drawn  attention  to  the  fact  that  most  images  are  not  smooth.  The 
PSI  method  is  predicated  on  locating  fe(x,  y)  in  the  correct  Lipschitz  space,  while  TV 
deblurring  assumes  fe(x,  y)  € BV(R2).  It  may  be  argued  that  such  refined  smoothness 
measures  are  primarily  applicable  to  f°°(x,y),  the  original  intensity  field  that  gave 
rise  to  the  digitized  finite  dimensional  object  fe(x,  y),  but  may  not  be  meaningful  for 
fe(x,y)  itself.  Indeed,  since  all  norms  are  equivalent  in  finite  dimensional  space,  it 
remains  to  be  seen  whether  such  abstruse  function  space  notions  are  ultimately  of  any 
computational  significance. 

Our  first  experiment  involves  a slightly  defocused  image  with  very  little  noise. 
The  original  sharp  USS  Eisenhower  image  is  shown  in  Figure  6(A).  Fourier  space 
multiplication  with  (27)  using  R = 0.03,  followed  by  the  addition  of  0.1%  multiplica- 
tive noise,  produced  the  blurred  Figure  6(B).  Because  of  the  low  noise  level,  we  chose 
P small  and  A large  in  (59),  as  recommended  in  [26].  With  /3  = 0.0001,  At  = O.l(Ax)2 
and  A = 300,  we  obtained  Figure  6(C)  at  T = 100 At.  Higher  values  of  A were  com- 
putationally unstable.  Moreover,  the  resulting  TV  image  did  not  improve  if  more 
time  steps  were  taken.  Figure  6(D)  is  the  PSI  deblurred  image  using  exact  values  for 
e,  M,  and  using  a = 0.362,  Cj  = 0.50,  from  Table  1.  Zooming  on  selected  parts  of 
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the  image  in  Figures  6(E)  and  6(F),  clearly  shows  significant  loss  of  structural  detail 
in  the  TV  image,  as  compared  with  PSI  deblurring.  For  completeness,  the  Ll  relative 
errors  in  this  experiment  were  as  follows:  true  Wiener  filtering  (not  shown)  1.67%  , 
PSI  method  2.16%,  TV  deblurring  6.83%,  and  Tikhonov-Miller  (not  shown)  4.71%. 

In  our  second  experiment,  the  sharp  English  Village  image  in  Figure  7(A)  was 
moderately  defocused  using  R = 0.06,  and  0.1%  multiplicative  noise  was  again  added 
to  form  Figure  7(B).  With  (3  and  A t as  in  Figure  6(C),  it  was  possible  to  choose 
A = 500,  and  obtain  Figure  7(C)  at  T — 100 At.  Again,  no  improvement  was  noted 
with  more  time  steps.  Figure  7(D)  is  the  PSI  deblurred  image  using  the  exact  values 
for  e,  i\/,  together  with  a = 0.439,  Cj  = 0.55,  from  Table  1.  Because  of  the  stronger 
blur,  more  information  is  now  lost  in  TV  deblurring.  Zooming  in  on  the  first  three 
houses  in  Figures  7(E)  and  7(F),  we  see  that  the  windows  and  roof  shingles  have 
been  almost  completely  eliminated  in  the  TV  image.  The  L 1 relative  errors  in  this 
experiment  were  as  follows:  true  Wiener  filtering  (not  shown)  1.98%,  PSI  method 
3.02%,  TV  deblurring  6.70%,  and  Tikhonov-Miller  (not  shown),  7.42%. 

In  Figure  4,  the  PSI  method’s  improvement  over  the  Tikhonov-Miller  method  can 
be  traced  to  the  fact  that  the  constraints  in  (29)  allowed  the  solution  to  be  too  rough. 
In  Figures  6 and  7,  PSI’s  improvement  over  the  total  variation  method  stems  from 
the  fact  that  the  minimum  principle  (57)  forces  the  solution  to  be  too  smooth.  Ap- 
parently, the  use  of  Lipschitz  spaces  to  calibrate  image  smoothness,  together  with  the 
direct  use  of  that  information  in  constraining  the  solutions  of  the  deblurring  problem, 
constitute  a significant  new  idea  in  image  deconvolution. 
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